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■ Abstract 

'-^ . We present a methodology to determine the best turbulence closure for 

an eddy-permitting ocean model: measurement of the error-landscape of the 
closure's subgrid spectral transfers and flux. Using a high-resolution bench- 
mark, we compare each closure's model of energy and enstrophy transfer 
to the actual transfer observed in the benchmark run. The error-landscape 
norms enable us to both make objective comparisons between the closures 
'. and to optimize each closure's free parameter for a fair comparison. We ap- 

^ ply this method to 6 different closures for forced-dissipative simulations of 

the barotropic vorticity equation on a f-plane (2D Navier-Stokes equation). 
The hyper-viscous closure most closely reproduces the enstrophy cascade es- 
pecially at larger scales due to the concentration of its dissipative effects 
to the very smallest scales. The viscous and Leith closures perform nearly 
as well especially at smaller scales where all three models were dissipative. 
The Smagorinsky closure dissipates enstrophy at the wrong scales. The an- 
ticipated potential vorticity closure was the only model to reproduce the 
inverse cascade, or back-scatter, of kinetic energy from the unresolved scales 
but would require high-order Laplacian corrections (to concentrate dissipa- 
tion only in the small scales) to be competitive. The Lagrangian-averaged 
a— model closure did not perform successfully for this fluid system: small- 
scale filamentation is preserved by the model but small-scale roll-up is pre- 
vented which inhibits the effects of diffusion. 

Keywords: 

Mesoscale eddies. Turbulent transfer. Parameterization, Oceanic 
turbulence, Eddy viscosity. Accuracy, Enstrophy 



Preprint submitted to Ocean Modeling 



July 26, 2012 



1. Introduction 



Turbulence closure models are required in the dynamical cores of global 
ocean-climate simulations. While grand challenge coupled climate simula- 
tions can use an ocean resolution of 0.1°(~ 10 km) to simulate timescales 
of decades, resolving the turbulent cascade for submesocale, 0(1 km), ed- 
dies remains computationally un achievable. For this reason, mesosca le ocean 
large-eddy simulations (MOLES; Fox-Kemper and Menemenlisl2008 ) are em- 
ployed. What is needed is a method to objectively compare the various 
proposed MOLES closures. Such a method is presented here, albeit in an 
idealized framework meant to serve as a basis for the evaluation and devel- 
opment of closures applicable to World Ocean simulations. 

Often, the closure approach taken is to set the dissipation scale equal to 
the grid scale. This is equivalent to setting the grid-scale Reynolds number 
to unity and is accomplished by simply using a constant viscosity, z/, that is 
much larger than the physical value (~ 10"^m^s"^) so that a well-resolved 
numerical simulation results. These large viscosities, however, also result in 
unphysical damping of the large scales. To reduce this effect while remaining 
in the paradigm of a linear dissipative model, the order of the Laplacian, 
A = V^, can be increased to = or higher. Such hyper- viscous mod- 
els are more scale-selective, applying dissipation concentrated near the grid 
scale (a new dissipation scale is derived from dimensional analysis of the A" 
dissipation and this scale is set equal to the grid scale). Turbulence is far 
more than a dissipative phenomenon, however, and purely dissipative mod- 
els cannot reproduce up-scale energy transfers due to interactions between 
scales. 

Another approach is to use what is known about turbulent cascades and 
apply dissipation only where it is re quired with a spati o-temporally vary- 
ing v iscosity, e.g., the Smagorinsky (jSmagorinskyl . Il963[ ) and Leith (ILeithl . 
19961 ) models. In the Smagorinsky model, the global average energy dissi- 



pation (due to a spatially uniform viscosity) is equated to the local dissi- 
pation at the grid scale because the turbulence is assumed homogeneous. 
The expression for i/*(a:,t) then follows from the 3D turbulence spectrum 
and dimensional analysis. For Leith, enstrophy dissipation and the 2D tur- 
bulence spectrum in an enstrophy cascade are used to derive the appro- 
priate i'^,(x,t) . How ever, the assumption of homogeneity is controversial 
( iLesieur et all . 120051) and there are also issues w i th vorticity d i ssipat ion at 
the boundaries (IFox-Kemper and Pedloskyl . l2004l : iFox-Kemperl . 120051 ). Yet, 
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the Leith model has been successful in improving; numer i cal sta bihty in global 
eddy-permitting models (IFox-Kemper and Menemenlisl . l2008l ). 

In 2D turbulent systems where enstrophy is clearly the quantity cas- 
cading to unresolved scales, methods to dissipate potential enstrophy while 
conserving energy have merit. This is accomplished by modifying the Cori- 
olis force in the momentum equation such that the transport of poten- 
tial vorticity is appropriately d iffusive while still being energetically neutral 
( ISadourny and Basdevantlll985l ). The anticipated potential vorticity method 
(APVM) reproduces both the physical tr ansfer of energy to lar ger scales and 
the dissipation of small-scale enstrophy (IVallis and Hual. Il988l). APVM has 
also been extended to variable- resolution grids ( IChen et al.l . l201ll ). However, 
it requires a high-order Laplacian correction to concentrate the eddy viscosity 
to the smallest scales ( IVallis and Hual . llQSSi ) . 

A recent approach is to use a mathematical regularization of the un- 
derlying equations, which ensures smooth (hence, computabl e) solutions 
a s th e closure model: e.g., the Lagrangian-averaged a— model (IHolm et al 



19981 ). It is dispersive rather than dissipative: the transport is by a spatially- 
smoothed velocity field (filter width ~ a). For three-dimensional (3D) in- 
compressible, non- rotating, and non-stratified flows the a— model does not 
prod uce sizeable computatiorial gain s because it develops rigid bodies in the 



flow (jPietarila Graham et al.l . 120071 1. This limitation can be overcome by 
applying the model to systems with a body force. It has been used suc- 
cessfuUy where there is a Lorentz f orce, in electrically conducting fluids 



(jPietarila Graham et al.l . l2009l 1201 ll ). and where there is a Coriolis force. 



in rotating fluids, e.g., the t wo-dimensional (2D) barotropic vorticity equa- 
tion (BVE) on a /3— plane (INadiga and Margolinl. 20011 : iHolm and Nadigal . 



20031 ). the shallow water equations flWingatd . 12004) . and a two-layer quasi- 



geostrophic (QG) model ( iHolm and Wingatd . l2005l ). 

For 2D flows, r elevant to this paper , the a — model enhances the inverse 
cascade of energy (INadiga and ShkoUerl . 1200 ll ) and in the enstrophy cascade 
regime, the rough kinetic energy and enstroph y spectra rema i n unc hanged 
{k~^ and k~^, respec tively) in the li r nit a — )■ oo (ILunasin et al.l . l2007l ). With 
forcing scaled to a, iLunasin et al.l ( l2007l ) found that increasing a lead to 
increasing the amount of fine structure and, consequently, to the need for 
increased resolution. They posited that with forcing unsealed, computational 
gains (instead of losses) might be realized. We will test whether or not this 
is so. 

The challenge in evaluating the effectiveness of LES closures for MOLES 
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should already be clear. Not only do many possible closures exist, but these 
closures often differ at the conceptual level of how unresolved turbulent mo- 
tion should be modeled. As such, we expect that the various possible closures 
will each excel in some plausible evaluation metric. The challenge is then to 
determine an approach, i.e an evaluation framework, that is both unbiased 
and fairly measures the effectiveness of the various closures in mimicking the 
influence of unresolved scales. The goal of this contribution is to do exactly 
that. 

Our approach here is to begin with the simplest system that we believe 
might be applicable to MOLES, with the understanding that the results ob- 
tained in such idealized systems will have to be reevaluated as the system 
complexity and realism increases. With this caveat in mind, we solve the 2D 
barotropic vorticity equation (2D BVE) in a doubly-periodic domain. The 
motivation for using the 2D BVE is to exploit the similarity of the the QG 
vorticity equation to the 2D BVE at spatial scales smaller than the Rossby ra- 
dius of deformati on. At these small scales, the QG vorticity equation reduces 
to the 2D BVE (jVallisl . l2006l ). The QG vorticity equation has a potential 
enstrophy cascade of QG eddies below the scale of the baroclinic instability. 
Similarly, the 2D BVE has an enstrophy cascade below the forcing scale, 
which serves here as an analog to the scale of the baroclinic instability. Fur- 
thermore, the robust analysis of spectral fluxes of energy and enstrophy in 
3D systems is sufficiently complicated to warrant starting at a lower spatial 
dimension. Since the 2D BVE system lacks the process of baroclinic insta- 
bility to initiate the turbulent mixing, we use large-scale, slowly varying in 
time, wind stress to activate the turbulence. As used in Ocean General Cir- 
culation Models (OGCMs), quadratic bottom-drag is used to obtain realistic 
equilibrium solutions. 

Details of the enstr ophy casca d e pro c ess can be measure d using enstro- 



phy transfer analysis (IKraichnanl . Il971t iMaltrud and Vallisl . 119931 ). Since 



a subgrid parameterization is envisioned to model all interactions between 
the modeled and missing scales, enstrophy transfer analysis provides a pow- 
erful metric to determine parameterization performance. Using the error- 
landscape of enstrophy flux, we quantify the performance of the six popular 
MOLES closures discussed above (the two linear dissipative models and the 
four nonlinear models derived from hypotheses about turbulence) employing 
a single, exponentially convergent, n umerical model, the Geophysical High 
Order Suite for Turbulence f GHOST: [Mininni et al.ll201lh . 

To compare the models, we start by computing a fully-resolved numerical 
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solution of a flow with a flxed, physical viscosity as the benchmark. This elim- 
inates the possibility of any bias between the parameterizations that could 
result from using any single MOLES at higher resolution as the benchmark. 
In Section [2l enstrophy transfer analysis is reviewed, and its application to 
MOLES and how this will be combined with the error-landscape is given. 
In Section [31 the details of the parameterizations are introduced and each 
parameterization is optimized with the error-landscape technique in order to 
make a fair comparison. 



2. Theory 

2.1. 2D turbulence 

For scales much smaller than the deformation radius, the qua si-geostroph ic 



potential vorticity equation reduces to the 2D-BVE (see, e.g., IVallisI l2006l ) 
The 2D-BVE are 



dtC + {tp,C} = F + uV^C - ■ V X (|n|u) 

u = -V X (#) , (1) 

where ( is the vorticity, tp the stream function, u the 2D velocity, F an 
external time-varying forcing to mimic wind stress, u the viscosity, z the 
out-of-plane unit vector, and Cn/h the coefficient of quadratic bottom drag. 
As a constant Coriolis parameter has no effect on 2D motion, Eqs. ([1]) also 
describe the 2D-BVE on a /—plane. Thus for scales much smaller than the 
deformation radius, where we intend to apply MOLES to ocean models, Eqs. 
dl]) are a very good approximation of the dynamics. 



A general overview of 2D turbulence theory (see, e.g., IVallisll2006[ ) is here 
presented (see Fig. [1]). Kinetic energy, ||up, and hence enstrophy, ^|CP, are 
injected into the fluid. Because both are quadratic ideal invariants (conserved 
in the absence of forcing and viscosity) and ^ = z ■ V x u, enstrophy cascades 
to smaller scales and energy undergoes an inverse cascade to larger scales 
(Fjortoft's theorem). Under the assumption of spectral locality, forcing and 
dissipation only affect the flow over a flnite range of scales and both cascades 
must have a constant flux (Fig. [H lower panel). The constant flux cascade 
regimes are called inertial ranges because only the inertial terms, u ■ Vu for 
energy and u ■ VC = {V^, C} for enstrophy, are non- negligible. Dimensional 
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inverse energy cascade 




Figure 1: Cartoon depicting 2D turbulence theory: kinetic energy spectrum (E(k), top 
panel) and fluxes (bottom panel) of enstropliy (n5(fc), blue solid line) and energy (n(fc), 
red dashed line). Kinetic energy undergoes an inverse cascade to large scales (negative 
flux) at the kinetic energy injection rate, e. Enstrophy undergoes a direct cascade to small 
scales at the enstrophy injection rate, 77. 
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analysis after equating constant fluxes to the inertial terms for energy and 
enstrophy yields a k~^/^ energy spectrum in the inverse energy cascade and 
a energy spectrum in the enstrophy cascade (Fig. [H upper panel). Fine 
theoretical detail s such as the logarithmic corr ection to the spectrum 
( iKraichnanl . Il97l[ ) and arguments about locality (IXiao et al.l . l2008l ) have here 
been omitted. 



2.2. Transfer analysis 

The {if), C} term is the only non-negligible term in the enstrophy cascade 
regime. It will also be shown in Section [273] to be the term whose small-scale 
interactions we need to parameterize. It is thus the focus of our comparison 
methodology. Other terms in the analysis will heretofore be abbreviated 
as J-" for forcing, T) for dissipation, and Q for large-scale drag (where, e.g., 
J-" = CP)- The time evolution of enstrophy at any physical position is given 
by the enstrophy-balance equation, 

dt\e = (dtC = -C{^, Cj + J' + V+Q. (2) 

The time evolution of the enstrophy spectrum at wavenumber k, Z{k), is 
similarly, 

dtZ{k) = CdtC = Sik) + T{k) + V{k) + Qik) , (3) 

where S{k) is the enstrophy transfer function (i.e., net enstrophy received by 
wavenumber k from all other wavenumbers) , 

S{k) = -C{k){M{k), (4) 

and where the Fourier transform is represented by " and complex-conjugation 
by ■*. The flux of enstrophy through wavenumber k, i.e., the sum of enstrophy 
given by all wavenumbers < k to wavenumbers > k (i.e., moving to smaller 
scales), is given by 

Us{k) = - S{k')dk' , (5) 

that is, the total amount of enstrophy flowing past wavenumber k to larger 
wavenumbers. The divergence of the flux is the transfer, S. 
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2.3. MOLES 

To reduce computational cost, MOLES solve only the largest scales of a 
flow. The remaining unresolved scales are filtered out. The filtering operation 
is indicated by ^ and the resulting equations are 



C 



dtC + {^A, C} = +'^V\ + a + F + vV'X - X (|n|u 



(6) 



where we have defined the subgrid term a = — {^jCI + IV^jC}- The subgrid 
term is the effects on the resolved scales by unresolved fluid motions. (Note 
that (T = V X V ■ r where r = — uu + uu, the traditional 3D LES subgrid 
stress tensor.) The time evolution of the enstrophy spectrum is now given 
by 

dtZ{k) = S{k) + L{k) + J^{k) + V{k) + Q{k) , (7) 

where S{k) is enstrophy received by wavenumber k from all other resolved 
wavenumbers, 

S{k) = -C{^X}. (8) 
and L{k) is the enstrophy received from all unresolved wavenumbers, 



m = ca. 



(9) 



The sum of enstrophy received from resolved and unresolved wavenumbers 
in a MOLES aims to be equal to the enstrophy transfer function from a fully 
resolved system, 

S{k) + L{k)^S{k), (10) 

for all wavenumbers smaller than the filter wavenumber. The flux of enstro- 
phy through wavenumber k due to resolved and modeled interactions is given 
by 



Urik) 



S{k') + L{k')dk' . 



2.4. Objective method: error-landscape of enstrophy flux 

To objectively cor npare parameterizations, we make use of th e error- 
landscape assessment ( iMeyers et al.l . l2003l . l2006l 120071 : iMeyersl . I2OIII ) in three 
bandwidths: large scale, intermediate scale and small scale. In a variation 
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on iMeyerg (120111 ). we employ a normalized Li norm error for these three 
wavelength bands expressed as, 



'° /;:;-;""in.(.)i* • 

to see how the model is performing in different bandwidths, where p = 1, 2, 3 
corresponds to large, intermediate and small scales, respectively. is the 
linear resolution and k = N{2 + 1)/100 is approximately the highest resolved 
wavenumber. The optimal parameter value for each method is in the region 
where all three of these error norms are near their minimum. (The term 
landscape is intuitive for two error norms.) Inter-model comparisons are 
then made using the Dp. 

2.5. Design of numerical experiments 

We employ a well-tested parallelized pseudo-spectral code ( iMininni et al.l . 
201ll ). The computational box has size [2-7r]^, and wave numbers vary from 



kmin = 1 to kmax = N/?) usiug a Standard 2/3 de-aliasing rule, where is the 
number of grid points per direction. To cast our results in meaningful units, 
the results are dimensionalized by / = /q/', t = tot' where indicates non- 
dimensionalized pseudo-spectral result and Iq = 504 x lO^/vr m and to = 1-2 x 
lO^s. To spin up our runs we begin with a 1008^ simulation (dimensionalized 
grid spacing 10 km) initialized with a few large-scale Fourier modes. The 
forcing is designed to mimic wind-stress at /c = 4: 

F = A{t) cos (4y + (py) - cos (4x + (p^) , (13) 

where (px = 7rsin(1.2 x 10~^s~^t) and (f)y = 7rsin(1.2 x lO^^vr s^^t/3) so 
that the wind varies with a period of about 60 days. The coefficient A is 
dynamically controlled to hold a stead enstrophy injection rate of 1.75 x 
10"^*^ s~^ to reduce the amount of required statistics to measure a constant 
flux cascade, i.e., 

= 1-75 X 10"^* . (14) 

JdA 

Time step is 600 s, u = 88m^s~^, and Co/h = 1.25 x 10~^m~^. The resulting 
root-mean-squared velocity is Vj-ms = 2.6 ms~^ and the forcing scale {k = 4) 
is Lp = 2520 km. The corresponding forcing-scale turnover time is 11 days 
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and the Reynolds number is Re = VrmsLp/j^ ~ 75, 000. This is run for over 
1300 days. The final turbulent state of this run is used as initial conditions 
for the benchmark and MOLES runs at u = 1.375 m^s~^. 



3. Analysis of parameterizations 

The goal of MOLES is to anticipate higher resolution results at an afford- 
able resolution by representing the effects of the unresolved eddies. To avoid 
any bias between the parameterizations, we use as the benchmark a fully 
resolved direct numerical solution (DNS) at a resolution of 8192^ of a flow 
with u = 1.375 m^s"^ Each MOLES is then run at a resolution of 1008^ and 
tested for its ability to reproduce the benchmark. This allows us to test the 
models' representations against a known solution: a DNS flow. Accordingly, 
the MOLES simulations also must use u = 1.375 m^s~^ in addition to the 
subgrid term or they should be compared, instead, to a z/ = benchmark 
which cannot be produced. 

The benchmark is run for 390 days, Vrms = 2.6 ms~^ and the corre- 
sponding forcing-scale turnover time is 11 days. The Reynolds number is 

4.8 X 10^. A snapshot of the vorticity of the benchmark run is shown 
in the Upper Left panel of Fig. [2l There are several large vortices of both 
signs. Over time, vortices stretch and fold vortex fllaments into the flne-scale 
features as seen. This is the enstrophy cascade process. This simulation is 
completely resolved and this cascade is arrested at the smallest scales by 
dissipation (Upper Right panel in Fig. [2]). Energy is injected by the forcing 
term (Lower Right panel in Fig. [2]) at a constant injection rate: an inverse 
cascade of energy and direct cascade of enstrophy result. The quadratic drag 
term serves to arrest the inverse cascade of kinetic energy and primarily re- 
moves energy (and enstrophy) at the largest scales. Though, it does remove 
both from a wide range of scales (Lower Left panel in Fig. [2]). 

The flux and resulting enstrophy spectrum for the benchmark are shown 
in Fig. [3l A power-law spectrum, Z{k) ~ is observed in the enstrophy 

cascade inertial range. It is steeper than the predicted spectrum due to 
the quadratic drag which acts at all scales of the flow: the difference between 
the enstrophy flux (solid line) and a constant flux is exactly the cumulative 
drag (dotted hne) . This st e eper spectrum is similar to the result for linear 
drag ( Danilov and Gurarid . 2001 ). Note that dissipation is not signiflcant 



for wavenumbers, k < 300. Reproducing this flow at a resolution of 1008^ 
{kmax = 336) will thus be a rigorous test of the parameterizations. 
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Figure 2: 8192^ benchmark; snapshot at 390 days for (Upper Left) vorticity, C, with 
thresholds ±1.5 x 10^^ s^^ (counter-clockwise vorticity is shown in yellow; clockwise in 
red); (Upper Right) absolute value of vorticity tendency due to dissipation, vS/^C,, black 
pixels arc 2.25 x 10~^s~^; (Lower Left) vorticity tendency due to quadratic drag, — x 
(|u|u), with thresholds ±1.38 x 10~^s~^; (Lower Right) vorticity tendency due to forcing, 
F, with thresholds ±1 x lO^'^s^^. 
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Figure 3: Benchmark run: (Top) Enstrophy flux (ns(fc), solid) and cumulative enstro- 
phy injection (dash-triple-dotted), dissipation (dashed), and quadratic drag (dotted). As 
quadratic drag operates at all but the dissipative scales, a constant enstrophy flux range is 
not seen. (Bottom) Compensated enstrophy spectrum, kZ{k), versus wavenumber, fc, for 
8192^ BVE benchmark. Quadratic drag acts at all scales and precludes a pure Zik) ~ 
spectrum. 
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3.1. Linear viscous parameterizations and their performance 

The simplest parameterization approximation is to assume the main ef- 
fect of subgrid turbulence is dissipative. Accordingly, the viscosity is often 
increased until a numerically resolved solution is possible. The subgrid term, 
cr, in the MOLES equation, Eq. (jS]) is then 

a={u'- , (15) 

with v' ^ V. A slightly more sophisticated approach is to add higher-order 
dissipation, hyper-viscosity, e.g. 

a = z/4V^C, (16) 

or even higher order. We focus on and parameterizations here. 

We apply the error-landscape of enstrophy flux technique to optimize the 
viscous model. The modeled flux, 11^ (A;), for the viscous model is shown in 
Fig. m Note that as the viscosity is varied, the modeled flux brackets both 
sides of the benchmark flux. This suggests an optimal v' for the model should 
be indicated by the enstrophy flux error-landscape. Indeed, the intermediate 
scale error, Di, has its minimum for u + v' = 11 m^s~^ where both Dq and D2 
are also near their minimal values. This is the optimal viscous model which 
we will compare to the other parameterizations. 

The approximate reproduction of the benchmark flux is accomplished by 
the action of the subgrid enstrophy transfer L{k) (Fig. \5^. As expected, the 
action of the viscous model is solely dissipative. The solid black line indicates 
what the true transfer with the unresolved scales should be, S{k) — S{k), 
calculated by filtering the benchmark down to a resolution of 1008^. The 
viscous model dissipates enstrophy over a much larger range of scales. More- 
over, since energy is dissipated as ~ v' Z{k) ~ /c~^'^, energy is unphysically 
dissipated at large scales. What the unresolved scales should be doing is 
contributing to the inverse cascade of energy as shown by the solid, black 
benchmark line. The enstrophy spectra are shown in Fig. |6l The result of 
too little dissipat i on is t he piling of small-scale thermal noise in the spectrum 



flCichowlas et all . 120051 ). 

By looking at the hyper-viscous model's flux error-landscape norms (Fig. 
[7]), we identify z/4 = 1.1 x 10^ m^s~^ as the optimal hyper- viscous model. The 
hyper-viscous model much more closely models the dissipation of enstrophy 
due to the unresolved scales than the viscous model, see Fig. |8l Addi- 
tionally, as the energy dissipation is ~ k'^Z ~ k^-^., energy is dissipated at 
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Figure 4: Viscous model: (Top) Modeled flux, nT(fc), for v + v' ^ 5.5m^s~^ (red 
dotted), llm^s"^ (green dashed), 16.5 m^s"^ (blue dash-dotted), 22m^s^^ (pink dash- 
triple-dotted), and 24.75 m^s~^ (cyan long-dashed) and IVs{k) for 8192^ BVE bench- 
mark (solid black). (Bottom) Flux error-landscape norms for large Dq (solid), inter- 
mediate Di (dashed), and small D2 (dot-dashed) resolved scales. The optimal value is 
V + v' = 11 m^s~^. 
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Figure 5: Viscous model: subgrid transfers for enstrophy (L(k\ Top) and energy {L{k)/k'^ , 
Bottom) and S{k) — S{k) for benchmark (solid black). The model is solely dissipative of 
enstrophy and energy. Exact viscosities arc denoted in Fig. HI 
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Figure 6: Viscous model: Compensated enstrophy spectrum; exact viscosities are denoted 
in Fig. m 



small scales (and at a much slower rate-note the difference in vertical scales 
for energy transfer in Figs. |5]and|H])- This is a marked improvement, but 
no solely-dissipative parameterization will model the inverse energy cascade 
mechanism. 

3.2. Leith model 



The Leith model is derived by dimensional analysis (jLeithl . Il996l ). The 
local enstrophy dissipation rate is estimated as 

?7* = z/*|V*C|\ (17) 

and an enstrophy cascade spectrum is assumed, 

Z{k) oc Tf'^k-^ . (18) 

The viscous range, k, is when the viscous enstrophy losses in a given wavenum- 
ber band, J vk"^ Z{k)dk, are comparable to the enstrophy injection, 77, or 

ri ~ u^k^ . (19) 

Setting the global average dissipation, i^, to the local, grid-scale dissipation 
rate, z/*, and equating Eqs. f|T7|) and f|T9|) . we find 

z/, oc |VC|(Ax)^ (20) 
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Figure 7: Hyper-viscous model: Flux error-landscape norms for large Dq (solid), inter- 
mediate Di (dashed), and small D2 (dot-dashed) resolved scales. The optimal value is 
= 1.1 X lO^m'^s^^ 



The BVE with Leith model, is ( iLeithl . Il996l : iFox-Kemper and Menemenlisl . 
2008h 

dtC + a = V ■ z/VC + V ■ v^VC , (21) 



where u 
is then 



for an infinite Reynolds number fiow. The Leith subgrid term 



V 



■j |VC|VC 



V TT 



(22) 



The subgrid transfers for the Leith model are very similar to the viscous 
model results (see Fig. |9]). This is to be expected as the Leith model is also 
solely-dissipative. What is surprising, however, is that there is strong enstro- 
phy dissipation at the forcing scale. This can be understood by looking at 
Fig. Uni The Leith viscosity oc |VC| and, therefore, is concentrated along 
the borders between oppositely-signed vortices. These large-scale coherent 
structures of enhanced dissipation then project on the small wavenumber 
Fourier-modes (bottom left panel of Fig. [TOj) . 

3.3. Smagorinsky model 



The Smagorinsky model ( ISmagorinskyl . Il963l : iLillyl . Il967( ) is the 3D pre- 
cursor of the Leith model. It is derived with a similar dimensional analysis 
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Figure 8: Hyper-viscous model; Subgrid transfers for enstropliy {L{k), Top) and en- 
ergy {L{k)/P, Bottom) for 1^4 = 2.2 x lO^m'^s-i (red dotted), 3.3 x lO^m'^s-i (green 
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Figure 9: Lcith model: Subgrid transfers for enstrophy {L{k), Top) and energy {L{k)/k'^, 
Bottom) for A ~ 0.5 (red dotted), A = 0.75 (green dashed), A = 1 (blue dash-dotted), 
A = 1.25 (pink dasli-triple-dotted), A = 1.5 (eyan long-dashed), and benchmark (black 
solid). The optimal model is A = 1. 



19 




Figure 10: Snapshots at 4 x lO^min for Lcith model, A = 1, i^, (Top Left) and Fourier 
power spectrum of (Bottom Left). Smagorinsky, As = 1, v^, (Top Right) and vorticity 
field (Bottom Right) are shown for reference. 
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as in Sec. 13. 2^ but assuming a 3D direct cascade of energy. Consequently, 
the model for eddy-viscosity is 



AAx\ 2 

(23) 



TT 



where Sij = {djVi + diVj)/2. For isotropi c, homogeneous 3D turbule nce the 



Smagorinsky Constant, C5 = A/7r^0.2 (iMeneveau and Katzl . l2000l ). 

The enstrophy flux and enstrophy spectrum for Smagorinsky (Fig. [TT]), 
highlight the fact that good spectra can be produced without necessarily 
reproducing the correct dynamics. The best spectra are produced for A5 = 
0.5 (blue dash-dotted) while the best flux is produced by A5 = 0.1 (red 
dotted). This is opposed to the case for the viscous model where the best flux 
and spectrum occur for u+u' = 11 m^s~^. The reason for the disparity is that 
the viscous parameterization captures the most important physical process, 
small-scale enstrophy dissipation, while the Smagorinsky model unphysically 
removes enstrophy (and energy) from the largest scales (see Fig. [12] and 
the real-space visualization of z/* in Fig. [TO]) . Therefore, even when the 
combination of modeling and numerical error produces a good spectrum, the 
model is not capturing the correct physical dynamics. 

3.4- Anticipated vorticity method (AVM) 



AVM (APVM when applied to potential vorticity, ISadourny and Basdevant 



19851 ) is so-called because it can be seen as substituting the forward-in-time 

-{^,Cn}, (24) 



vorticity in the BVE, 

Cn+l ~ Cn 



9 

where 9 is the time step for the anticipation. Substituting this anticipated 
value, Cn+i in Eq- dl]) results in the lowest-order AVM, 

5tC = -{^,Cn} + ^M{^,Cj}. (25) 

In practice, to weight the subgrid model to smaller scales, 

^ = -T^MV^-{^,C}}, (26) 

max 

In this study we have used m = 1 as even this order of diffusive opera- 
tor is not practical in finite-volume and finite-difference schemes typically 
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Figure 11: Smagorinsky model: (Top) Modeled flux, Urik), for As = 0.1 (red dotted), 
As = 0.3 (green dashed), A5 = 0.5 (blue dash-dotted), A5 = 1 (pink dash-triple-dotted), 
and A5 = 2 (cyan long-dashed) and lls{k) for 8192^ BVE benchmark (soUd black). (Bot- 
tom) Compensated enstrophy spectrum. 
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Figure 12: Smagorinsky model: subgrid transfers for enstrophy, L{k). The model dis- 
sipates enstrophy (and energy) unphysically from the large scales. Exact viscosities are 
denoted in Fig. [TTl 



used in global ocean modeling. AVM is not Galilean invariant, i.e., it does 
not conserve momentum, but it exactly conserves energy while dissipating 
enstrophy. 

As AVM dissipates enstrophy at small scales, L{k) < for large k (see 
Fig. [T3|) . it must also remove some small-scale energy, k~'^L{k) < 0. Since 
AVM exactly conserves energy, this energy shows up at large scales. AVM 
is the only parameterization studied here that reproduces this signature of 
the correct transfer. The physical effect, however, is over estimated by at 
least an order of magnitude. This can be mitigated by reducing 9. However, 
too small 9 (0.1 2 5dt fo r our flow) results in an excess of energy at all scales 
(jVallis and Hual . Il988l ). For m = 1, as used here, AVM is unable to mimic 
that eddy viscosity should only act in a small range of wavenumbers near 
kmax (jVallis and Hual . Il988l ). For the AVM, the optimal value of 9 is difficult 
to choose. Fig. [HI For the smallest resolved scales, 9 = 0.25, is best, but the 
largest resolved scales are best reproduced by 6* = 1. We compromise and 
choose 9 = 0.5. 
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Figure 13: AVM: Subgrid transfers of cnstrophy {L{k), Top) and energy {L{k)/k'^, Bot- 
tom) for 6* = (NO MODEL, red dotted), ^ 0.15625 (green dashed), 0.25 (blue dash- 
dotted), 0.5 (pink dash-triple-dotted), and 1 (cyan long-dashed). The subgrid model 
transfer in AVM changes sign so that the model dissipates no energy, sum of i(fc)/fc^ over 
all wavenumbers is o(10^^^), while cnstrophy dissipation (sum of L{k)) is o(l). The neg- 
ative energy dissipation at large scales mimics the inverse cascade from unresolved scales, 
though too strongly. 
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Figure 14: AVM: (Top) Modeled flux, IVrik). Exact values of 9 are given in Fig. [131 
(Bottom) Flux error-landscape norms. Optimal value is chosen to be = 0.5. 
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3.5. a— model 

The a— model takes a different approach than the other parameteriza- 
tions. It is a non-dissipative, solely dispersive model - a mathematical regu- 
larization (smooth, and hence compu table solutions are ensured even in the 



limit z/ — )■ 0) of the fluid equations (IHolm et al.l . Il998l ). The result is that 
the vorticity is advected by a smoothed velocity, ui' = (1 — a^V^)~^u, with 
a (second) filter scale ~ a, 

dtC + V-{w,c)=^^^C + :F + V+Q, (27) 

where V ■ (^"^Cj — {^s? C}- The alpha subgrid term is 

(x = MC}-{^.,C}- (28) 

The subgrid transfers. Fig. [151 for the a— model are very large and in 
the wrong direction. As the model dissipates neither energy nor enstrophy 
the transfers are conservative; they remove energy and enstrophy from above 
the forcing scale and deposit them below the forcing scale. As the filter 
width, a, is increased so is the amount of large-scale energy and enstrophy 
moved down-sc ale to scales larger t han a (vertical lines in Fig. [T5|l . This is 



consistent with lLunasin et al.l ( 120071 ) 's finding that more numerical resolution 
is required as a is increased. 

The physical affect of the a— model is visualized in Fig. [161 small-scale 
vortical motions are removed from the advecting field. The central feature, 
the yellow(light) V-shaped counter-clockwise vorticity rolls up less as a is 
increased. The filament continues to be stretched to smaller scales, however, 
by the larger-scale motions. The combination of these two effects means 
that there are still many thin filaments (small-scales), but without the small- 
scale winding they are less well mixed. This reduces the effect of dissipation 
and hence the increased resolution requirements found in this and the earlier 
study. 

3.6. Comparison of parameterizations 

The subgrid transfers of the six parameterizations are compared in Fig. 
[T7i The hyper-viscous is closest to mimicking the true transfer, though for 
the largest wavenumbers, k > 200, the viscous and Leith parameterizations 
perform similarly. The AVM is the only method that reproduces the correct 
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Figure 15: a— model: subgrid transfers of enstrophy {L{k), Top) and energy {L{k)/k'^, 
Bottom) for a — Ax (red dotted). 2 Ax (green dashed), 9 Ax (blue dash-dotted; vertical 
line shows wavenumber), 16 Ax (pink dash-triple-dotted; vertical line shows wavenumber), 
and benchmark (solid black). Due to numerical cancellation noise in Eq. (|28| . smoothing 
has been applied to the plots. 
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Figure 16: Vorticity field, (: Vortex merger event (tracked to center of field of view 
which is 1/5^ of the entire domain). Time runs from top to bottom starting Omin after 
initialization in steps of min. 1st column LANS a — 2 Ax, 2nd column 9 Ax, 3rd 
column 16Aa;. a = 2Ax is the most realistic result. 
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sign of the energy transfer, but it removes enstrophy preferentially from scales 
around 6Ax instead of the smallest resolved scale s. This method would likely 
perform better for m > 1 (jVallis and Hual . Il988[ ). 

The enstrophy flux error landscape norms are given in Fig. [181 The 
a— model and Smagorinsky are the obvious outliers with poor performance 
at all scales. At the largest resolved scales, the hyper-viscous model has 
the best performance due to its extremely limited effect at these scales. At 
intermediate resolved scales, 100 < k < 200, the viscous and hyper-viscous 
parameterizations have very similar performance. For the smallest resolved 
scales, the viscous method is the best. At all scales, the Leith model has 
performance very close to the viscous model. The AVM also has similar per- 
formance except at the smallest scales. Again, this could likely be mitigated 
by using a larger value of m. 

This similarity in performance between hyper-viscous, viscous, and Leith 
parameterizations can also been seen in the resulting enstrophy spectra. Fig. 
[T9l All three methods have the same spectra for k ^ 100. The a— model 
increases the need for resolution, opposite of what is needed for a turbulence 
parameterization, and so does not reduce the pile-up of thermalization noise 
in the small scales. As seen in the previous results, the AVM method with 
m = 1 is dissipative at too large scales to perform as well as the viscous, 
hyper-viscous, or Leith parameterizations. Smagorinsky performs poorly be- 
cause it removes enstrophy from the largest rather than the smallest resolved 
scales. 



4. Conclusions 

We have compared six popular turbulence parameterizations in the en- 
strophy cascade regime of the barotropic vorticity equation on a /—plane 
(equivalently, 2D Navier-Stokes) in forced-dissipative simulations. The hyper- 
viscous, viscous, and Leith models all perform well down to about 10 Ax. The 
hyper- viscous model reproduces the largest-resolved-scales {1 < k < 100) flux 
the best of the three and the viscous model best reproduces the smallest- 
resolved-scales {k > 200) flux. The Leith model is expected to carry-over 
its performance to anisotropic flows (e.g., the 3D baroclinic ocean system) 
which would be challenging for the viscous and hyper-viscous models. The 
Smagorinsky model does not work in the enstrophy cascade regime-it re- 
moves enstrophy from the largest rather than the smallest resolved scales. 
The anticipated vorticity method without a strong enough weighting to small 
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Figure 17: Subgrid transfers of enstrophy {L{k), Top) and energy {L{k)/k'^, Bottom) 
for benchmark (solid black), hyper-viscous V4 = 1.1 x 10~^m'^s~^ (red dotted), LANS 
a = 2A.T (green dashed), viscous i' = llm^s"^ (blue dash-dotted), Leith A = 1 (pink 
dash-triple-dotted); AVM 6 = 0.5 (cyan long-dashed), and Smagorinsky A5 =0.1 (solid 
grey). 
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Figure 18: Enstrophy flux error-landscape norms for hyper-viscous {1^4), a— model (a), 
viscosity (v), Leith (A), AVM (6), and Smagorinsky (A5): large Dq (solid), intermediate 
Di (dashed), and small D2 (dot-dashed) resolved scales. 




Figure f9: Enstrophy spectra for benchmark (solid black), hyper-viscous 1^4 = 1.1 x 
10~^m^s~^ (red dotted), LANS a = 2Ax (green dashed), viscous 1/ — llm^s"^ (blue 
dash-dotted), Leith A = 1 (pink dash-triple-dotted), AVM 6 = 0.5 (cyan long-dashed), 
and Smagorinsky As — 0.1 (solid grey). 
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scales, higher-orders than acting on the anticipated vorticity, does not 
perform as well as the prior three parameterizations. As even this order of 
diffusive operator is not practical in the finite-volume and finite-differenc e 
schemes typically used in global ocean modeling (e.g., iRingler et al.ll2010l ). 
we chose not investiga te higher-order s. 

We have confirmed iLunasin et al.l (120071 ) 's suggestion that the Lagrangian- 
averaged a— model does not perform as a turbulence model in this system. 
The reason is two-fold. First, large-scale motions continue to stretch vortex 
filaments creating more and more small scales. Second, small-scale vorticity 
is smoothed away from the advecting velocity field by the model and fila- 
ments no longer roll up. The result is small-scale features that are harder to 
dissipate: the opposite of what one looks for in a turbulence model. 

One possible MOLES closure has not been scoped here, the use of mono- 
tone transport as the model for LES closure. These closures, commonly 
referred to Monotone Implicit Large-Eddy Simulation (MILES), require the 
evaluation of the nonlinear transport be carried out in physical space, some- 
thing that is not possible within the global spectral model utilized for this 
study. Our future work, discussed briefly below, will utilize a traditional 
finite-volume approach where an evaluation of MILES will be possible. 



Subgrid transfers have been measured before, e.g., for the APVM (jVallis and Hua , 



19881 ). an d again for the APVM , hyper- viscosity, and implicit large-eddy sim- 



ulations (IThuburn et al.l . l201ll ). Error-lan dscapes for LES have been pro- 



duced for var ious quantities like spectra (IMeyers et al.1 . l2003l . l2006l . 12007 



Meyersl . l201ll ). By combining these two techniques, however, we have in- 



troduced a method for determining the optimal turbulence parameterization 
also in flows different than those considered here: the error-landscape of the 
enstrophy flux at small-scales in a 2D flow can be replaced by the error- 
landscape of the modeled potential enstrophy flux in a 3D baroclinic system. 
We emphasize that MOLES comparisons based on spectra alone do not en- 
sure that the correct dynamics are being reproduced by a parameterization. 
For example, consider the As = 0.5 {Cs ~ 0.16) result for the Smagorinsky 
model (blue dash-dotted line in Fig. [TT]). The spectrum is best approximated 
by this run, but for the wrong reasons as the non-linear flux is more poorly 
reproduced than for As = 0.1. For the viscous model, however, which phys- 
ically correctly removes enstrophy from the small scales, both the spectrum 



and the flux are best reproduced for v 



11m s . In this latter case, 



the spectrum is reproduced because the dynamics are reproduced. 

For a 3D baroclinic system, at the scales on which the MOLES acts, the 
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system will be approximately 2D. We then expect that our results will hold: 
Smagorinsky, the a— model, and APVM with m = 1 will not perform as well 
as viscosity, hyper- viscosity, and Leith. In fact, because of the anisotropic 
diffusion offered by Leith, it will likely perform the best. Our results will not 
extend to the 3D system, however, if additional physics comes in to play, like 
vertical mixing over small horizontal scales. 

Our next step is to move into an idealized, 3D baroclinic system, most 
likely a re-entrant zonal channel that can serve as an idealized Antarctic 
Circumpolar Current. While the move to three dimensions allows for the 
direct simulation of baroclinic instability, it also necessitates the develop- 
ment of analysis tools that can accurately account for energy and enstrophy 
transfers between the disparate horizontal and vertical scales of motion. Fur- 
thermore, the move to a 3D baroclinic system will also likely entail the use 
of a height-based vertical coordinate. Such a system requires the transport 
of one or more tracers in order to close the system via an equation of state. 
As a result, the evaluation of turbulent mixing is no longer confined to the 
momentum equation, but now also must include the mixing of tracers. And 
finally, as we move to more realistic and, thus, bounded domains, our ability 
to simulate the governing equations, as well as analyze the fluxes, via global 
spectral expansions will be increasingly cumbersome. As a result, we plan 
on utilizing a traditional finite- volume global ocean model in the next phase 
of this study. We will also l i kely e xplore the idea of evaluating fluxes based 
on a local expansion (lEyinkl . l2005l ) instead of a global expansion. 
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